* This file: .\replication\figures_2_3_4_5_6\figures_2_3_4_5_6.do

* Run this file first, to generate data-sets:
* ./replication/figures_2_3_4_5_6/data_for_figures_2_3_4_5_6.do
* Need to dowload data-files first: 
* Freddie Mac’s "Single Family Loan-Level Dataset" (http://www.freddiemac.com/research/datasets/sf_loanlevel_dataset.page, last accessed June 3, 2021 [https://perma.cc/KU32-9TXH]), Q1-1999 to Q4-2019

**********
** Figure 2: Gini at year-month level, full sample
**********
use ./mortgages_usa/year_month_full.dta, clear
keep if year<=2019

foreach dv in gini cv { 
	reg `dv' c.myear
	cap drop `dv'_poly1
	predict `dv'_poly1 if e(sample)
	label var `dv'_poly1 "Linear fit-line"
	reg `dv' c.myear##c.myear
	cap drop `dv'_poly2
	predict `dv'_poly2 if e(sample)
	label var `dv'_poly2 "Quadratic fit-line"
	reg `dv' c.myear##c.myear##c.myear
	cap drop `dv'_poly3
	predict `dv'_poly3 if e(sample)
	label var `dv'_poly3 "Cubic fit-line"
	reg `dv' c.myear##c.myear##c.myear##c.myear
	cap drop `dv'_poly4
	predict `dv'_poly4 if e(sample)
	label var `dv'_poly4 "Fourth degree polynomial fit-line"
}

foreach dv in gini cv {
	tw  ///
		(scatter `dv'  myear, sort ms(o) mc(gs7)) ///
		(line `dv'_poly1 myear, sort lc(black) lp(dash)) ///
		(line `dv'_poly3 myear, sort lc(black) lp(solid)) ///
		, xline(576 577 578 579, lc(gs12) lw(medthick) lp(solid)) xlabel(#5) scheme(538bw) legend(off) ytitle("`: var label `dv''") xlabel(#20, angle(90) labsize(vsmall)) ysize(4) xsize(4) xtitle("") 	note("", size(vsmall)) ylabel(, labsize(vsmall))
	graph save   figure_2_`dv'.gph
	graph export figure_2_`dv'.pdf
	graph export figure_2_`dv'.svg
}


**********
** Figure 3: FICO scores in sample, over time
**********
use ./mortgages_usa/year_month_full.dta, clear
keep if year<=2019

cap drop _merge
merge 1:1 myear using ./mortgages_usa/year_month_full_fico620.dta
drop _merge
keep if year<=2019

gen fico_range=fico_max-fico_min
label var fico_range "Range of FICO-scores"
label var fico_mean "Mean of FICO-scores"
label var fico_p50 "Median of FICO-scores"

keep myear fico_mean fico_p50 fico_range fico620
rename fico620 fico_620
reshape long fico_, i(myear) j(q) string

gen quantity=""
replace quantity="Percentage with FICO<620" if q=="620"
replace quantity="Mean FICO score" if q=="mean"
replace quantity="Median FICO score" if q=="p50"
replace quantity="Range of FICO scores" if q=="range"

tw  ///
	(scatter fico_  myear, sort ms(o) mc(gs7)) ///
	, xline(576 577 578 579, lc(gs12) lw(medthick) lp(solid)) xlabel(#5) scheme(538bw) legend(off) ytitle("") xlabel(#20, angle(90) labsize(vsmall)) ysize(4) xsize(4) xtitle("") 	note("", size(vsmall)) ylabel(, labsize(vsmall)) ///
	by(quantity, r(2) yrescale note(""))

graph save   figure_3.gph
graph export figure_3.pdf
graph export figure_3.svg


******************************************
* Analysis at the year-month-FICO2d level
******************************************
* Figure 4: Interest rate spread over time (year-month-FICO-2d level)
* Figure 5: Interest rate spread over time (year-month-FICO2d level), at FICO2d
* Figure 6: RD estimate at FICO-2d-levels, different polynomials

use ./mortgages_usa/year_month_sample_fico2d.dta, clear

sum N, detail
drop if N<100 // groups are too small
drop if fico_2d>=82 // groups are too small
sum N
di `r(N)'*`r(mean)'

foreach v in gini cv sd my dev1 dev2 {
	replace `v'=. if `v'==0
}

gen post=(myear>=576)
foreach dv in gini cv dev1 dev2  {
	reg `dv' c.myear##i.post##i.fico_2d
	cap drop `dv'_poly1
	predict `dv'_poly1 if e(sample)
	label var `dv'_poly1 "Linear fit-line"
	reg `dv' c.myear##c.myear##i.post##i.fico_2d
	cap drop `dv'_poly2
	predict `dv'_poly2 if e(sample)
	label var `dv'_poly2 "Quadratic fit-line"
	reg `dv' c.myear##c.myear##c.myear##i.post##i.fico_2d
	cap drop `dv'_poly3
	predict `dv'_poly3 if e(sample)
	label var `dv'_poly3 "Cubic fit-line"
	reg `dv' c.myear##c.myear##c.myear##c.myear##i.post##i.fico_2d
	cap drop `dv'_poly4
	predict `dv'_poly4 if e(sample)
	label var `dv'_poly4 "Fourth degree polynomial fit-line"
}

foreach dv in gini cv dev1 dev2  {
	reg `dv' c.myear##i.post
	cap drop `dv'_apoly1
	predict `dv'_apoly1 if e(sample)
	label var `dv'_apoly1 "Linear fit-line"
	reg `dv' c.myear##c.myear##i.post
	cap drop `dv'_apoly2
	predict `dv'_apoly2 if e(sample)
	label var `dv'_apoly2 "Quadratic fit-line"
	reg `dv' c.myear##c.myear##c.myear##i.post
	cap drop `dv'_apoly3
	predict `dv'_apoly3 if e(sample)
	label var `dv'_apoly3 "Cubic fit-line"
	reg `dv' c.myear##c.myear##c.myear##c.myear##i.post
	cap drop `dv'_apoly4
	predict `dv'_apoly4 if e(sample)
	label var `dv'_apoly4 "Fourth degree polynomial fit-line"
}

***********************************************************************
* Figure 4: Interest rate spread over time (year-month-FICO-2d level)
***********************************************************************

foreach dv in gini cv dev1 dev2 {
	tw  ///
		(scatter `dv'  myear, jitter(1) sort ms(o) msize(tiny) mc(gs7)) ///
		(line `dv'_apoly1 myear if myear< 576, sort lw(medthick) lc(black) lp(solid)) ///
		(line `dv'_apoly1 myear if myear>=576, sort lw(medthick) lc(black) lp(solid)) ///
		(line `dv'_apoly2 myear if myear< 576, sort lw(medthick) lc(black)   lp(shortdash)) ///
		(line `dv'_apoly2 myear if myear>=576, sort lw(medthick) lc(black)   lp(shortdash)) ///
		(line `dv'_apoly3 myear if myear< 576, sort lw(thick) lc(gs4) lp(dash)) ///
		(line `dv'_apoly3 myear if myear>=576, sort lw(thick) lc(gs4) lp(dash)) ///
		, xline(576 577 578 579, lc(gs12) lw(medthick) lp(solid)) xlabel(#5) ///
		scheme(538bw) ///
		ytitle("`: var label `dv''") xlabel(#20, angle(90) labsize(vsmall)) ysize(4) xsize(4) xtitle("") 	note("", size(vsmall)) ylabel(, labsize(vsmall)) ///
		legend(order(2 4 6) r(1) pos(6))
	graph save   figure_4_`dv'.gph
	graph export figure_4_`dv'.pdf
	graph export figure_4_`dv'.svg
}

***********************************************************************
* Figure 5: Interest rate spread over time (year-month-FICO2d level), at FICO2d
***********************************************************************

foreach dv in gini cv dev1 dev2 {
	tw  ///
		(scatter `dv'  myear, sort ms(o) msize(tiny) mc(gs7)) ///
		(line `dv'_poly1 myear if myear< 576, sort lw(medthick) lc(black) lp(solid)) ///
		(line `dv'_poly1 myear if myear>=576, sort lw(medthick) lc(black) lp(solid)) ///
		(line `dv'_poly2 myear if myear< 576, sort lw(medthick) lc(black)   lp(shortdash)) ///
		(line `dv'_poly2 myear if myear>=576, sort lw(medthick) lc(black)   lp(shortdash)) ///
		(line `dv'_poly3 myear if myear< 576, sort lw(thick) lc(gs4) lp(dash)) ///
		(line `dv'_poly3 myear if myear>=576, sort lw(thick) lc(gs4) lp(dash)) ///
		, xline(576 577 578 579, lc(gs12) lw(medthick) lp(solid)) xlabel(#5) ///
		scheme(538bw) ///
		ytitle("`: var label `dv''") xlabel(#20, angle(90) labsize(vsmall)) xtitle("") 	note("", size(vsmall)) ylabel(, labsize(vsmall)) ///
		legend(order(2 4 6) r(1) pos(6)) ///
		by(fico_2d, c(5) note("")) ///
		ysize(6) xsize(5)

	graph save   figure_5_`dv'.gph
	graph export figure_5_`dv'.pdf
	graph export figure_5_`dv'.svg
}


* net install rdrobust, from(https://raw.githubusercontent.com/rdpackages/rdrobust/master/stata) replace
* net install rdlocrand, from(https://raw.githubusercontent.com/rdpackages/rdlocrand/master/stata) replace
* net install rddensity, from(https://raw.githubusercontent.com/rdpackages/rddensity/master/stata) replace
* net install lpdensity, from(https://raw.githubusercontent.com/nppackages/lpdensity/master/stata) replace
* net install rdmulti, from(https://raw.githubusercontent.com/rdpackages/rdmulti/master/stata) replace
* net install rdpower, from(https://raw.githubusercontent.com/rdpackages/rdpower/master/stata) replace


***********************************************************************
* Figure 6: RD estimate at FICO-2d-levels, different polynomials
***********************************************************************
* Also: all RD estimates
clear all

cap log close
log using ./tables_rd/tables_rd_estimates.txt, text

di "** Current data **"
di "`= c(current_date)'"
about
which rdrobust

use ./mortgages_usa/year_month_sample_fico2d.dta, clear
drop if N<100
drop if fico_2d>=82 // groups are too small
foreach v in gini cv sd my dev1 dev2 {
	replace `v'=. if `v'==0
}
gen post=(myear>=576)

// potential covariates
levelsof month, local(M)
foreach m of local M {
	gen month_`m'=(month==`m')
}
levelsof year, local(Y)
foreach y of local Y {
	gen year_`y'=(year==`y')
}
recode month (1/6=1) (7/12=2), gen(half)
cap drop timehalf
egen timehalf=group(year half), label
numlabel timehalf, add
tab timehalf

recode month (1/3=1) (4/6=2) (7/9=3) (10/12=4), gen(quarter)
recode month (1/2=1) (3/4=2) (5/6=3) (7/8=4) (9/10=5) (11/12=6), gen(bimonth)

xtset fico_2d myear
foreach dv in gini cv dev1 dev2  {
	cap drop l1_`dv'
	gen l1_`dv'=l1.`dv'
}


cd ./tables_rd/

estimates clear
foreach dv in gini {
	foreach p in 1 2 3 {
		forvalues j=62(1)81 {
			eststo s0_`p'_`dv'_`j': rdrobust `dv'  myear if fico_2d==`j', c(576) p(`p')
			eststo c1_`p'_`dv'_`j': rdrobust `dv'  myear if fico_2d==`j', c(576) p(`p') covs(month_1 month_2 month_3 month_4 month_5 month_6 month_7 month_8 month_9 month_10 month_11)
			eststo c2_`p'_`dv'_`j': rdrobust `dv'  myear if fico_2d==`j', c(576) p(`p') covs(l1_`dv')
			eststo c3_`p'_`dv'_`j': rdrobust `dv'  myear if fico_2d==`j', c(576) p(`p') covs(l1_`dv' month_1 month_2 month_3 month_4 month_5 month_6 month_7 month_8 month_9 month_10 month_11)
		}
	}
}

local t s0
foreach p in 1 2 3 {
	esttab s0_1_gini_?? using figure_6_`t'_`p'_gini.rtf ///
		,  ///
		stats(ci_rb kernel N p q N_l N_r N_h_l N_h_r h_l h_r b_l b_r, fmt(0 0 0 0 0 0 0 0 0 3 3) label("Robust 95% CI" "Kernel Type" "Observations" "Order Loc. Poly. (p)" "Order Bias (q)" "Number of obs (left of cutoff)" "Number of obs (right of cutoff)" "Eff. Number of obs (left of cutoff)" "Eff. Number of obs (right of cutoff)" "BW est. (left of cutoff)" "BW est. (right of cutoff)" "BW bias (left of cutoff)" "BW bias (right of cutoff)")) ///
		addnotes("at FICO-2d levels; s0: no covariates; c1: months; c2: lagged DV; c3: months + lagged DV") 
}

cd ..

* Figure 6: RD estimate at FICO-2d-levels, different polynomials
coefplot s0_1_gini_??, pstyle(p1) bylabel(Linear fit) || s0_2_gini_??, bylabel(Quadratic fit) || s0_3_gini_??, bylabel(Cubic fit) scheme(538bw) nokey  byopts(r(1)) ///
	ylabel( ///
	0.545 "620-629" ///
	0.593 "630-639" ///
	0.641 "640-649" ///
	0.689 "650-659" ///
	0.737 "660-669" ///
	0.785 "670-679" ///
	0.833 "680-689" ///
	0.881 "690-699" ///
	0.929 "700-709" ///
	0.977 "710-719" ///
	1.025 "720-729" ///
	1.073 "730-739" ///
	1.121 "740-749" ///
	1.169 "750-759" ///
	1.217 "760-769" ///
	1.265 "770-779" ///
	1.313 "780-789" ///
	1.361 "790-799" ///
	1.409 "800-809" ///
	1.457 "810-819" ///
	) ///
	ytitle("FICO score band") ///
	xtitle("Estimated coefficient with 95% CIs") ///
	xlabel(0(.5)2.5)
	
	graph save   figure_6.gph
	graph export figure_6.pdf
	graph export figure_6.svg


cd ./tables_rd/

* Table 1: Regression discontinuity estimates
	foreach p in 1 2 3 4 {
		eststo pn_p`p': rdrobust gini myear, c(576) p(`p') all kernel(tri) bwselect(mserd) 
		eststo pc_p`p': rdrobust gini myear, c(576) p(`p') all kernel(tri) bwselect(mserd) ///
			vce(cluster fico_2d)
	}
	esttab pn_p? pc_p?, star mtitle ///
		stats(ci_rb kernel N p q N_l N_r N_h_l N_h_r h_l h_r b_l b_r, fmt(0 0 0 0 0 0 0 0 0 3 3) label("Robust 95% CI" "Kernel Type" "Observations" "Order Loc. Poly. (p)" "Order Bias (q)" "Number of obs (left of cutoff)" "Number of obs (right of cutoff)" "Eff. Number of obs (left of cutoff)" "Eff. Number of obs (right of cutoff)" "BW est. (left of cutoff)" "BW est. (right of cutoff)" "BW bias (left of cutoff)" "BW bias (right of cutoff)")) 
	esttab pc_p* using table_1.rtf,  ///
		star mtitle se onecell ///
		varwidth(10) modelwidth(5) ///
		stats(ci_l_cl ci_r_cl ci_l_rb ci_r_rb ci_rb bwselect kernel p q N N_l N_r N_h_l N_h_r h_l h_r b_l b_r cmd depvar runningvar vce_select clustvar , fmt(2 2 2 2 0 0 0 0 0 0 0 0 0 0 2 2) ///
		label("95% CI, LB" "95% CI, UB" "Robust 95% CI, LB" "Robust 95% CI, UB" "Robust 95% CI" "BW type" "Kernel" "Order Loc. Poly. (p)" "Order bias (q)" "N" "N (l)" "N (r)" "Eff. N (l)" "Eff. N (r)" "BW est. (l)" "BW est. (r)" "BW bias (l)" "BW bias (r)" "Command" "Outcomevar" "Runningvar" "VCE method" "Clustervar" )) ///
		addnotes("`e(cmd)' `e(depvar)' `e(runningvar)', c(`e(c)') p(`e(p)') all kernel(`e(kernel)') bwselect(`e(bwselect)') vce(`e(vce_select)' `e(clustvar)')")

* Table D6: Different bandwidth selection procedures
	* Different bandwidth selection
	foreach m in ///
		mserd /// one common MSE-optimal bandwidth selector for the RD treatment effect estimator.
		msetwo /// two different MSE-optimal bandwidth selectors (below and above the cutoff) for the RD treatment             effect estimator.
		msesum  /// one common MSE-optimal bandwidth selector for the sum of regression estimates (as opposed to             difference thereof).
		msecomb1  /// for min(mserd,msesum).
		msecomb2  /// for median(msetwo,mserd,msesum), for each side of the cutoff separately.
		cerrd  /// one common CER-optimal bandwidth selector for the RD treatment effect estimator.
		certwo  /// two different CER-optimal bandwidth selectors (below and above the cutoff) for the RD treatment 		effect estimator.
		cersum  /// one common CER-optimal bandwidth selector for the sum of regression estimates (as opposed to 		difference thereof).
		cercomb1  ///  for min(cerrd,cersum).
		cercomb2  ///  for median(certwo,cerrd,cersum), for each side of the cutoff separately.
		{
		eststo bn_`m': rdrobust gini myear, c(576) p(1) all kernel(tri) bwselect(`m') 
		eststo bc_`m': rdrobust gini myear, c(576) p(1) all kernel(tri) bwselect(`m') ///
			vce(cluster fico_2d)
	}
	esttab bn_*, star ///
		stats(ci_rb kernel N p q N_l N_r N_h_l N_h_r h_l h_r b_l b_r, fmt(0 0 0 0 0 0 0 0 0 3 3) label("Robust 95% CI" "Kernel Type" "Observations" "Order Loc. Poly. (p)" "Order Bias (q)" "Number of obs (left of cutoff)" "Number of obs (right of cutoff)" "Eff. Number of obs (left of cutoff)" "Eff. Number of obs (right of cutoff)" "BW est. (left of cutoff)" "BW est. (right of cutoff)" "BW bias (left of cutoff)" "BW bias (right of cutoff)")) compress
	esttab bc_*, star ///
		stats(ci_rb kernel N p q N_l N_r N_h_l N_h_r h_l h_r b_l b_r, fmt(0 0 0 0 0 0 0 0 0 3 3) label("Robust 95% CI" "Kernel Type" "Observations" "Order Loc. Poly. (p)" "Order Bias (q)" "Number of obs (left of cutoff)" "Number of obs (right of cutoff)" "Eff. Number of obs (left of cutoff)" "Eff. Number of obs (right of cutoff)" "BW est. (left of cutoff)" "BW est. (right of cutoff)" "BW bias (left of cutoff)" "BW bias (right of cutoff)")) compress

	esttab bc_* using table_d6_bandwidth.rtf,  ///
		star mtitle se onecell ///
		varwidth(10) modelwidth(5) ///
		stats(ci_l_cl ci_r_cl ci_l_rb ci_r_rb ci_rb bwselect kernel p q N N_l N_r N_h_l N_h_r h_l h_r b_l b_r cmd depvar runningvar vce_select clustvar , fmt(2 2 2 2 0 0 0 0 0 0 0 0 0 0 2 2) ///
		label("95% CI, LB" "95% CI, UB" "Robust 95% CI, LB" "Robust 95% CI, UB" "Robust 95% CI" "BW type" "Kernel" "Order Loc. Poly. (p)" "Order bias (q)" "N" "N (l)" "N (r)" "Eff. N (l)" "Eff. N (r)" "BW est. (l)" "BW est. (r)" "BW bias (l)" "BW bias (r)" "Command" "Outcomevar" "Runningvar" "VCE method" "Clustervar" )) ///
		addnotes("`e(cmd)' `e(depvar)' `e(runningvar)', c(`e(c)') p(`e(p)') all kernel(`e(kernel)') bwselect(`e(bwselect)') vce(`e(vce_select)' `e(clustvar)')")

* Table D7: Different kernel functions
	* Different kernel functions
	foreach m in ///
		triangular ///
		epanechnikov ///
		uniform {
		eststo kn_`m': rdrobust gini myear, c(576) p(1) all kernel(`m') bwselect(mserd) 
		eststo kc_`m': rdrobust gini myear, c(576) p(1) all kernel(`m') bwselect(mserd) ///
			vce(cluster fico_2d)
	}
	esttab kn_*, star mtitle se ///
		stats(ci_rb kernel N p q N_l N_r N_h_l N_h_r h_l h_r b_l b_r, fmt(0 0 0 0 0 0 0 0 0 3 3) label("Robust 95% CI" "Kernel Type" "Observations" "Order Loc. Poly. (p)" "Order Bias (q)" "Number of obs (left of cutoff)" "Number of obs (right of cutoff)" "Eff. Number of obs (left of cutoff)" "Eff. Number of obs (right of cutoff)" "BW est. (left of cutoff)" "BW est. (right of cutoff)" "BW bias (left of cutoff)" "BW bias (right of cutoff)")) compress
	esttab kc_*, star mtitle se  ///
		stats(ci_rb kernel N p q N_l N_r N_h_l N_h_r h_l h_r b_l b_r, fmt(0 0 0 0 0 0 0 0 0 3 3) label("Robust 95% CI" "Kernel Type" "Observations" "Order Loc. Poly. (p)" "Order Bias (q)" "Number of obs (left of cutoff)" "Number of obs (right of cutoff)" "Eff. Number of obs (left of cutoff)" "Eff. Number of obs (right of cutoff)" "BW est. (left of cutoff)" "BW est. (right of cutoff)" "BW bias (left of cutoff)" "BW bias (right of cutoff)")) compress

	esttab kc_* using table_d7_kernel.rtf,  ///
		star mtitle se onecell ///
		varwidth(10) modelwidth(5) ///
		stats(ci_l_cl ci_r_cl ci_l_rb ci_r_rb ci_rb bwselect kernel p q N N_l N_r N_h_l N_h_r h_l h_r b_l b_r cmd depvar runningvar vce_select clustvar , fmt(2 2 2 2 0 0 0 0 0 0 0 0 0 0 2 2) ///
		label("95% CI, LB" "95% CI, UB" "Robust 95% CI, LB" "Robust 95% CI, UB" "Robust 95% CI" "BW type" "Kernel" "Order Loc. Poly. (p)" "Order bias (q)" "N" "N (l)" "N (r)" "Eff. N (l)" "Eff. N (r)" "BW est. (l)" "BW est. (r)" "BW bias (l)" "BW bias (r)" "Command" "Outcomevar" "Runningvar" "VCE method" "Clustervar" )) ///
		addnotes("`e(cmd)' `e(depvar)' `e(runningvar)', c(`e(c)') p(`e(p)') all kernel(`e(kernel)') bwselect(`e(bwselect)') vce(`e(vce_select)' `e(clustvar)')")

* Table D8: Covariate adjusted estimates
	* Covariates
	estimates clear
	eststo cov_none: rdrobust gini myear, c(576) p(1) all kernel(tri) bwselect(mserd) ///
		vce(cluster fico_2d)
	eststo cov_months: rdrobust gini myear, c(576) p(1) all kernel(tri) bwselect(mserd) ///
		covs(month_*) vce(cluster fico_2d)
	eststo cov_lag: rdrobust gini myear, c(576) p(1) all kernel(tri) bwselect(mserd) ///
		covs(l1_gini) vce(cluster fico_2d) 
	eststo cov_int: rdrobust gini myear, c(576) p(1) all kernel(tri) bwselect(mserd) ///
		covs(my)  vce(cluster fico_2d)
	eststo cov_all: rdrobust gini myear, c(576) p(1) kernel(tri) bwselect(mserd) ///
		covs(month_* my l1_gini) vce(cluster fico_2d)

	esttab cov_*, star se mtitles ///
		stats(ci_rb kernel N p q N_l N_r N_h_l N_h_r h_l h_r b_l b_r, fmt(0 0 0 0 0 0 0 0 0 3 3) label("Robust 95% CI" "Kernel Type" "Observations" "Order Loc. Poly. (p)" "Order Bias (q)" "Number of obs (left of cutoff)" "Number of obs (right of cutoff)" "Eff. Number of obs (left of cutoff)" "Eff. Number of obs (right of cutoff)" "BW est. (left of cutoff)" "BW est. (right of cutoff)" "BW bias (left of cutoff)" "BW bias (right of cutoff)")) 

	esttab cov_* using table_d8_covariates.rtf,  ///
		star mtitle se onecell ///
		varwidth(10) modelwidth(5) ///
		stats(ci_l_cl ci_r_cl ci_l_rb ci_r_rb ci_rb bwselect kernel p q N N_l N_r N_h_l N_h_r h_l h_r b_l b_r cmd depvar runningvar vce_select clustvar , fmt(2 2 2 2 0 0 0 0 0 0 0 0 0 0 2 2) ///
		label("95% CI, LB" "95% CI, UB" "Robust 95% CI, LB" "Robust 95% CI, UB" "Robust 95% CI" "BW type" "Kernel" "Order Loc. Poly. (p)" "Order bias (q)" "N" "N (l)" "N (r)" "Eff. N (l)" "Eff. N (r)" "BW est. (l)" "BW est. (r)" "BW bias (l)" "BW bias (r)" "Command" "Outcomevar" "Runningvar" "VCE method" "Clustervar" )) ///
		addnotes("`e(cmd)' `e(depvar)' `e(runningvar)', c(`e(c)') p(`e(p)') all kernel(`e(kernel)') bwselect(`e(bwselect)') vce(`e(vce_select)' `e(clustvar)')")

* Table D12: Placebo outcomes
	* Different outcome variables
	estimates clear
	foreach dv of varlist gini ltv my fico_mean { // fico_p50 fico_mean fico_max   purpose_* channel_* proptype_* borrowers cltv dti ltv {
		eststo dv_`dv': rdrobust `dv' myear, c(576) p(1) all kernel(tri) bwselect(mserd) ///
			vce(cluster fico_2d)
	}
	esttab dv_*, star mtitle se ///
		stats(ci_rb kernel N p q N_l N_r N_h_l N_h_r h_l h_r b_l b_r, fmt(0 0 0 0 0 0 0 0 0 3 3) label("Robust 95% CI" "Kernel Type" "Observations" "Order Loc. Poly. (p)" "Order Bias (q)" "Number of obs (left of cutoff)" "Number of obs (right of cutoff)" "Eff. Number of obs (left of cutoff)" "Eff. Number of obs (right of cutoff)" "BW est. (left of cutoff)" "BW est. (right of cutoff)" "BW bias (left of cutoff)" "BW bias (right of cutoff)")) 

	esttab dv_* using table_d12_ouctomes.rtf,  ///
		star mtitle se onecell ///
		varwidth(10) modelwidth(5) ///
		stats(ci_l_cl ci_r_cl ci_l_rb ci_r_rb ci_rb bwselect kernel p q N N_l N_r N_h_l N_h_r h_l h_r b_l b_r cmd depvar runningvar vce_select clustvar , fmt(2 2 2 2 0 0 0 0 0 0 0 0 0 0 2 2) ///
		label("95% CI, LB" "95% CI, UB" "Robust 95% CI, LB" "Robust 95% CI, UB" "Robust 95% CI" "BW type" "Kernel" "Order Loc. Poly. (p)" "Order bias (q)" "N" "N (l)" "N (r)" "Eff. N (l)" "Eff. N (r)" "BW est. (l)" "BW est. (r)" "BW bias (l)" "BW bias (r)" "Command" "Outcomevar" "Runningvar" "VCE method" "Clustervar" )) ///
		addnotes("`e(cmd)' `e(depvar)' `e(runningvar)', c(`e(c)') p(`e(p)') all kernel(`e(kernel)') bwselect(`e(bwselect)') vce(`e(vce_select)' `e(clustvar)')")

* Table D11: Sensitivity to observations near the cutoff (donut hole approach)
	* Donut
	gen t0=myear-576
	tab myear if t0==0
	estimates clear
	forvalues j=0(1)3 {
		cap noisily eststo dn_`j': rdrobust gini t0 if abs(t0)>=`j', c(0) p(1) all kernel(tri) 
		cap noisily eststo dc_`j': rdrobust gini t0 if abs(t0)>=`j', c(0) p(1) all kernel(tri) ///
			vce(cluster fico_2d)
	}
	esttab dn_*, star mtitle se ///
		stats(ci_rb kernel N p q N_l N_r N_h_l N_h_r h_l h_r b_l b_r, fmt(0 0 0 0 0 0 0 0 0 3 3) label("Robust 95% CI" "Kernel Type" "Observations" "Order Loc. Poly. (p)" "Order Bias (q)" "Number of obs (left of cutoff)" "Number of obs (right of cutoff)" "Eff. Number of obs (left of cutoff)" "Eff. Number of obs (right of cutoff)" "BW est. (left of cutoff)" "BW est. (right of cutoff)" "BW bias (left of cutoff)" "BW bias (right of cutoff)")) compress
	esttab dc_*, star mtitle se ///
		stats(ci_rb kernel N p q N_l N_r N_h_l N_h_r h_l h_r b_l b_r, fmt(0 0 0 0 0 0 0 0 0 3 3) label("Robust 95% CI" "Kernel Type" "Observations" "Order Loc. Poly. (p)" "Order Bias (q)" "Number of obs (left of cutoff)" "Number of obs (right of cutoff)" "Eff. Number of obs (left of cutoff)" "Eff. Number of obs (right of cutoff)" "BW est. (left of cutoff)" "BW est. (right of cutoff)" "BW bias (left of cutoff)" "BW bias (right of cutoff)")) compress

	esttab dc_* using table_d11_donut.rtf,  ///
		star mtitle se onecell ///
		varwidth(10) modelwidth(5) ///
		stats(ci_l_cl ci_r_cl ci_l_rb ci_r_rb ci_rb bwselect kernel p q N N_l N_r N_h_l N_h_r h_l h_r b_l b_r cmd depvar runningvar vce_select clustvar , fmt(2 2 2 2 0 0 0 0 0 0 0 0 0 0 2 2) ///
		label("95% CI, LB" "95% CI, UB" "Robust 95% CI, LB" "Robust 95% CI, UB" "Robust 95% CI" "BW type" "Kernel" "Order Loc. Poly. (p)" "Order bias (q)" "N" "N (l)" "N (r)" "Eff. N (l)" "Eff. N (r)" "BW est. (l)" "BW est. (r)" "BW bias (l)" "BW bias (r)" "Command" "Outcomevar" "Runningvar" "VCE method" "Clustervar" )) ///
		addnotes("`e(cmd)' `e(depvar)' `e(runningvar)', c(`e(c)') p(`e(p)') all kernel(`e(kernel)') bwselect(`e(bwselect)') vce(`e(vce_select)' `e(clustvar)')")

* Table D13: Placebo cutoffs
	* Different/placebo cutoffs	
	tab myear if t0==0
	*estimates clear
	*forvalues j=6(1)18 {
	forvalues j=6(6)30 {
		cap noisily eststo Lcn_`j': rdrobust gini t0 if t0<0, c(-`j') p(1) all kernel(tri) 
		cap noisily eststo Lcc_`j': rdrobust gini t0 if t0<0, c(-`j') p(1) all kernel(tri) ///
			vce(cluster fico_2d)
		cap noisily eststo Rcn_`j': rdrobust gini t0 if t0>=0, c(`j') p(1) all kernel(tri) 
		cap noisily eststo Rcc_`j': rdrobust gini t0 if t0>=0, c(`j') p(1) all kernel(tri) ///
			vce(cluster fico_2d)
	}
	esttab Lcn_* Lcc_*, star mtitle ///
		stats(ci_rb kernel N p q N_l N_r N_h_l N_h_r h_l h_r b_l b_r, fmt(0 0 0 0 0 0 0 0 0 3 3) label("Robust 95% CI" "Kernel Type" "Observations" "Order Loc. Poly. (p)" "Order Bias (q)" "Number of obs (left of cutoff)" "Number of obs (right of cutoff)" "Eff. Number of obs (left of cutoff)" "Eff. Number of obs (right of cutoff)" "BW est. (left of cutoff)" "BW est. (right of cutoff)" "BW bias (left of cutoff)" "BW bias (right of cutoff)")) compress
	esttab Rcn_* Rcc_*, star mtitle ///
		stats(ci_rb kernel N p q N_l N_r N_h_l N_h_r h_l h_r b_l b_r, fmt(0 0 0 0 0 0 0 0 0 3 3) label("Robust 95% CI" "Kernel Type" "Observations" "Order Loc. Poly. (p)" "Order Bias (q)" "Number of obs (left of cutoff)" "Number of obs (right of cutoff)" "Eff. Number of obs (left of cutoff)" "Eff. Number of obs (right of cutoff)" "BW est. (left of cutoff)" "BW est. (right of cutoff)" "BW bias (left of cutoff)" "BW bias (right of cutoff)")) compress

	coefplot Lcn_* Lcc_*, scheme(538bw)
	coefplot Rcn_* Rcc_*, scheme(538bw)

	esttab dc_0 Lcc_* using table_d13_placebocutoff.rtf,  ///
		star mtitle se onecell ///
		varwidth(10) modelwidth(5) ///
		stats(ci_l_cl ci_r_cl ci_l_rb ci_r_rb ci_rb bwselect kernel p q N N_l N_r N_h_l N_h_r h_l h_r b_l b_r cmd depvar runningvar vce_select clustvar , fmt(2 2 2 2 0 0 0 0 0 0 0 0 0 0 2 2) ///
		label("95% CI, LB" "95% CI, UB" "Robust 95% CI, LB" "Robust 95% CI, UB" "Robust 95% CI" "BW type" "Kernel" "Order Loc. Poly. (p)" "Order bias (q)" "N" "N (l)" "N (r)" "Eff. N (l)" "Eff. N (r)" "BW est. (l)" "BW est. (r)" "BW bias (l)" "BW bias (r)" "Command" "Outcomevar" "Runningvar" "VCE method" "Clustervar" )) ///
		addnotes("`e(cmd)' `e(depvar)' `e(runningvar)', c(`e(c)') p(`e(p)') all kernel(`e(kernel)') bwselect(`e(bwselect)') vce(`e(vce_select)' `e(clustvar)')")

* Table D14: Masspoints
	* Masspoints
	foreach m in ///
		adjust ///
		off {
		eststo mn_`m': rdrobust gini myear, masspoints(`m') c(576) p(1) all kernel(tri) bwselect(mserd) 
		eststo mc_`m': rdrobust gini myear, masspoints(`m') c(576) p(1) all kernel(tri) bwselect(mserd) ///
			vce(cluster fico_2d)
	}
	esttab mn_*, star mtitle se ///
		stats(ci_rb kernel N p q N_l N_r N_h_l N_h_r h_l h_r b_l b_r, fmt(0 0 0 0 0 0 0 0 0 3 3) label("Robust 95% CI" "Kernel Type" "Observations" "Order Loc. Poly. (p)" "Order Bias (q)" "Number of obs (left of cutoff)" "Number of obs (right of cutoff)" "Eff. Number of obs (left of cutoff)" "Eff. Number of obs (right of cutoff)" "BW est. (left of cutoff)" "BW est. (right of cutoff)" "BW bias (left of cutoff)" "BW bias (right of cutoff)")) compress
	esttab mc_*, star mtitle se  ///
		stats(ci_rb kernel N p q N_l N_r N_h_l N_h_r h_l h_r b_l b_r, fmt(0 0 0 0 0 0 0 0 0 3 3) label("Robust 95% CI" "Kernel Type" "Observations" "Order Loc. Poly. (p)" "Order Bias (q)" "Number of obs (left of cutoff)" "Number of obs (right of cutoff)" "Eff. Number of obs (left of cutoff)" "Eff. Number of obs (right of cutoff)" "BW est. (left of cutoff)" "BW est. (right of cutoff)" "BW bias (left of cutoff)" "BW bias (right of cutoff)")) compress

	esttab mc_* using table_d14_masspoints.rtf,  ///
		star mtitle se onecell ///
		varwidth(10) modelwidth(5) ///
		stats(ci_l_cl ci_r_cl ci_l_rb ci_r_rb ci_rb bwselect kernel p q N N_l N_r N_h_l N_h_r h_l h_r b_l b_r cmd depvar runningvar vce_select clustvar , fmt(2 2 2 2 0 0 0 0 0 0 0 0 0 0 2 2) ///
		label("95% CI, LB" "95% CI, UB" "Robust 95% CI, LB" "Robust 95% CI, UB" "Robust 95% CI" "BW type" "Kernel" "Order Loc. Poly. (p)" "Order bias (q)" "N" "N (l)" "N (r)" "Eff. N (l)" "Eff. N (r)" "BW est. (l)" "BW est. (r)" "BW bias (l)" "BW bias (r)" "Command" "Outcomevar" "Runningvar" "VCE method" "Clustervar" )) ///
		addnotes("`e(cmd)' `e(depvar)' `e(runningvar)', c(`e(c)') p(`e(p)') all kernel(`e(kernel)') bwselect(`e(bwselect)') vce(`e(vce_select)' `e(clustvar)')")

* Table D9: Different specification of the running variable I
	* Level of analysis I (based on year-month-FICO-2d data)
	cap drop quarter
	recode month (1/3=1) (4/6=2) (7/9=3) (10/12=4), gen(quarter)

	cap drop trimester
	recode month (1/4=1) (5/8=2) (9/12=3), gen(trimester)

	cap drop biannual
	recode month (1/6=1) (7/12=2), gen(biannual)

	cap drop annual
	recode month (1/12=1), gen(annual)

	estimates clear
	foreach v in month quarter trimester biannual annual {
		cap drop t0
		cap drop time
		egen time=group(year `v'), label(time, replace)
		sum time if year==2008 & `v'==1
		gen t0=time-`r(mean)'
		di "`v'"
		tab time if t0==0

		cap noisily eststo ln_`v': rdrobust gini t0, c(0) p(1) all kernel(tri) 
		cap noisily eststo lc_`v': rdrobust gini t0, c(0) p(1) all kernel(tri) ///
				vce(cluster fico_2d)
	}
	esttab ln_*, star mtitle ///
		stats(ci_rb kernel N p q N_l N_r N_h_l N_h_r h_l h_r b_l b_r, fmt(0 0 0 0 0 0 0 0 0 3 3) label("Robust 95% CI" "Kernel Type" "Observations" "Order Loc. Poly. (p)" "Order Bias (q)" "Number of obs (left of cutoff)" "Number of obs (right of cutoff)" "Eff. Number of obs (left of cutoff)" "Eff. Number of obs (right of cutoff)" "BW est. (left of cutoff)" "BW est. (right of cutoff)" "BW bias (left of cutoff)" "BW bias (right of cutoff)")) compress
	esttab lc_*, star mtitle ///
		stats(ci_rb kernel N p q N_l N_r N_h_l N_h_r h_l h_r b_l b_r, fmt(0 0 0 0 0 0 0 0 0 3 3) label("Robust 95% CI" "Kernel Type" "Observations" "Order Loc. Poly. (p)" "Order Bias (q)" "Number of obs (left of cutoff)" "Number of obs (right of cutoff)" "Eff. Number of obs (left of cutoff)" "Eff. Number of obs (right of cutoff)" "BW est. (left of cutoff)" "BW est. (right of cutoff)" "BW bias (left of cutoff)" "BW bias (right of cutoff)")) compress


	esttab lc_* using table_d9_unit1.rtf,  ///
		star mtitle se onecell ///
		varwidth(10) modelwidth(5) ///
		stats(ci_l_cl ci_r_cl ci_l_rb ci_r_rb ci_rb bwselect kernel p q N N_l N_r N_h_l N_h_r h_l h_r b_l b_r cmd depvar runningvar vce_select clustvar , fmt(2 2 2 2 0 0 0 0 0 0 0 0 0 0 2 2) ///
		label("95% CI, LB" "95% CI, UB" "Robust 95% CI, LB" "Robust 95% CI, UB" "Robust 95% CI" "BW type" "Kernel" "Order Loc. Poly. (p)" "Order bias (q)" "N" "N (l)" "N (r)" "Eff. N (l)" "Eff. N (r)" "BW est. (l)" "BW est. (r)" "BW bias (l)" "BW bias (r)" "Command" "Outcomevar" "Runningvar" "VCE method" "Clustervar" )) ///
		addnotes("`e(cmd)' `e(depvar)' `e(runningvar)', c(`e(c)') p(`e(p)') all kernel(`e(kernel)') bwselect(`e(bwselect)') vce(`e(vce_select)' `e(clustvar)')")

* Table D10: Different specification of the running variable II
	* Level of analysis II (based on year-X-FICO-2d data)
	clear all
	estimates clear
	cd ..
	cd ./mortgages_usa/
	foreach v in month quarter trimester biannual annual {
		use year_`v'_sample_fico2d.dta, clear
		sum N
		drop if N<100
		gen tt="`v'"
		drop if year>2019
		
		tab `v'
		cap drop t0
		cap drop time
		cap label drop time
		egen time=group(year `v'), label(time, replace)
			sum time if year==2008 & `v'==1
		gen t0=time-`r(mean)'
		di "`v'"
		tab time if t0==0

		cap noisily eststo ln2_`v': rdrobust gini t0, c(0) p(1) all kernel(tri) 
		cap noisily eststo lc2_`v': rdrobust gini t0, c(0) p(1) all kernel(tri) ///
				vce(cluster fico_2d)
	}
	esttab ln2_*, star mtitle ///
		stats(ci_rb kernel N p q N_l N_r N_h_l N_h_r h_l h_r b_l b_r, fmt(0 0 0 0 0 0 0 0 0 3 3) label("Robust 95% CI" "Kernel Type" "Observations" "Order Loc. Poly. (p)" "Order Bias (q)" "Number of obs (left of cutoff)" "Number of obs (right of cutoff)" "Eff. Number of obs (left of cutoff)" "Eff. Number of obs (right of cutoff)" "BW est. (left of cutoff)" "BW est. (right of cutoff)" "BW bias (left of cutoff)" "BW bias (right of cutoff)")) compress
	esttab lc2_*, star mtitle ///
		stats(ci_rb kernel N p q N_l N_r N_h_l N_h_r h_l h_r b_l b_r, fmt(0 0 0 0 0 0 0 0 0 3 3) label("Robust 95% CI" "Kernel Type" "Observations" "Order Loc. Poly. (p)" "Order Bias (q)" "Number of obs (left of cutoff)" "Number of obs (right of cutoff)" "Eff. Number of obs (left of cutoff)" "Eff. Number of obs (right of cutoff)" "BW est. (left of cutoff)" "BW est. (right of cutoff)" "BW bias (left of cutoff)" "BW bias (right of cutoff)")) compress

	cd ..
	cd ./tables_rd/
	esttab lc2_* ln2_trimester using table_d10_unit2.rtf,  ///
		star mtitle se onecell ///
		varwidth(10) modelwidth(5) ///
		stats(ci_l_cl ci_r_cl ci_l_rb ci_r_rb ci_rb bwselect kernel p q N N_l N_r N_h_l N_h_r h_l h_r b_l b_r cmd depvar runningvar vce_select clustvar , fmt(2 2 2 2 0 0 0 0 0 0 0 0 0 0 2 2) ///
		label("95% CI, LB" "95% CI, UB" "Robust 95% CI, LB" "Robust 95% CI, UB" "Robust 95% CI" "BW type" "Kernel" "Order Loc. Poly. (p)" "Order bias (q)" "N" "N (l)" "N (r)" "Eff. N (l)" "Eff. N (r)" "BW est. (l)" "BW est. (r)" "BW bias (l)" "BW bias (r)" "Command" "Outcomevar" "Runningvar" "VCE method" "Clustervar" )) ///
		addnotes("`e(cmd)' `e(depvar)' `e(runningvar)', c(`e(c)') p(`e(p)') all kernel(`e(kernel)') bwselect(`e(bwselect)') vce(`e(vce_select)' `e(clustvar)')")

cap log close
cd ..
